//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

destring bwg , replace 

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark all sample
reg $baseline_child
gen sample0 = e(sample)

// mark baseline sample
reg $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen sample1 = e(sample)

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std) if sample1 == 1, method(adf);
predict parenting_age1 if sample1 == 1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std) if sample1 == 1, 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3 if sample1 == 1, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

global inputs_std   cum_avg_daycare_36m_sum parenting_ages13
keep   ihdp site tg bwg sample* $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs ppvtstd wasifsiq

// merge in ages 5 and 8
merge 1:1 ihdp using iqchildhoodihdp
keep  if _merge == 3
drop _merge

// impute iq
replace kidppvt5 = wippsif5 if kidppvt5 ==.
replace wippsif5 = kidppvt5 if wippsif5 ==.

replace ppvtstd8 = fsiq8    if ppvtstd8 ==.
replace fsiq8    = ppvtstd8 if fsiq8    ==.

replace ppvtstd  = wasifsiq if ppvtstd  ==.
replace wasifsiq = ppvtstd  if wasifsiq ==.

reg fsiq8 ppvtstd8 kidppvt5 wippsif5 wasifsiq ppvtstd
gen sample2 = e(sample) 

// outcomes
merge 1:1 ihdp using adultihdp
keep if _merge == 3
drop    _merge

// construct "positive" outcomes
replace idle18y         = 1 - idle18y
replace sch_eversped18y = 1 - sch_eversped18y
replace sch_math18y     = 1 - sch_math18y
replace sch_read18y     = 1 - sch_read18y
replace sch_thrp18y     = 1 - sch_thrp18y
replace teen18y         = 1 - teen18y

gen     smoke18 = .
replace smoke18 = 1 if cigs18y >  0  & cigs18y != .
replace smoke18 = 0 if cigs18y == 0 & cigs18y != .
replace smoke18 = 1 - smoke18

gen     absence18 = .
replace absence18 = 0 if sch_dabs18y < 5  & sch_dabs18y != .
replace absence18 = 1 if sch_dabs18y >= 5 & sch_dabs18y != .
replace absence18 = 1 - absence18 

// marking missing
global ed18  sch_eversped18y sch_math18y sch_read18y sch_test18y
global beh18 smoke18 idle18y sch_thrp18y teen18y
egen missing     = rowmiss($ed18 $beh18)

// average
egen avg_outcome_educ18 = rowmean($ed18)        if missing == 0
egen avg_outcome_beh18  = rowmean($beh18)       if missing == 0
egen avg_outcome18      = rowmean(avg_outcome_educ18 avg_outcome_beh18) if missing == 0

// define sample
reg     avg_outcome18 if sample2 == 1
replace sample2 = 0   if e(sample) == 0
replace sample2 = 0   if   sample1 == 0
keep if sample2 == 1

// strata
global perry_strata
global abc_strata
global ihdp_strata strata(bwg site treat)

rename avg_outcome_educ18 educ18
rename avg_outcome_beh18  beh18 
rename avg_outcome18      all18
rename sex male

foreach var of varlist educ18 beh18 all18 {
	reg     `var'   iqcage stndscor ppvtstd8 kidppvt5 wippsif5 fsiq8 ppvtstd wasifsiq
	predict `var'r, resid
	summ    `var'r if tg == 0
	replace `var'r = (`var'r - r(mean))/r(sd)
}

global sample1 if twin != . & male !=.
global sample2 if twin == 0 & male !=.
global sample3 if twin == 1 & male !=.
global sample4 if twin != . & male == 1
global sample5 if twin == 0 & male == 1
global sample6 if twin == 1 & male == 1
global sample7 if twin != . & male == 0 
global sample8 if twin == 0 & male == 0 
global sample9 if twin == 1 & male == 0 

// inference
foreach var of varlist educ18 beh18 all18 educ18r beh18r all18r {
	// all
	foreach sample of numlist 1(1)9 {
		gen `var'_s`sample' = `var' ${sample`sample'}
		bootstrap, strata(bwg site tg) reps($bootstraps): reg `var'_s`sample' tg ${sample`sample'}
		matrix b = e(b)
		matrix V = e(V)
		matrix se`var'_`sample' = sqrt(V[1,1])
		matrix md`var'_`sample' = b[1,1]
		matrix  c`var'_`sample' = b[1,2]
		matrix  t`var'_`sample' = abs(md`var'_`sample'[1,1]/se`var'_`sample'[1,1])
		matrix df`var'_`sample' = e(N) - e(rank)
		matrix  p`var'_`sample' = 2*(1 - normal(t`var'_`sample'[1,1]))
		matrix  n`var'_`sample' = e(N)
		matrix           `var'_`sample' = [c`var'_`sample',md`var'_`sample',p`var'_`sample',n`var'_`sample']
		matrix  colnames `var'_`sample' = c`var'_`sample' md`var'_`sample' p`var'_`sample' n`var'_`sample'
	}
	// twins vs. no twins
	foreach sample of numlist 2 3 5 6 8 9 {
		summ    `var'_s`sample'                             ${sample`sample'}  & tg == 0 
		replace `var'_s`sample' = `var'_s`sample' - r(mean) ${sample`sample'}
	}	
	egen `var'pool   = rowtotal(`var'_s2 `var'_s3), m 
	egen `var'male   = rowtotal(`var'_s5 `var'_s6), m
	egen `var'female = rowtotal(`var'_s8 `var'_s9), m
	
	foreach sample in pool male female { 
		bootstrap, strata(bwg site tg) reps($bootstraps): reg `var'`sample' twin if tg == 1
		matrix b = e(b)
		matrix V = e(V)
		matrix se`var'_`sample' = sqrt(V[1,1])
		matrix md`var'_`sample' = b[1,1]
		matrix  t`var'_`sample' = abs(md`var'_`sample'[1,1]/se`var'_`sample'[1,1])
		matrix df`var'_`sample' = e(N) - e(rank)
		matrix  p`var'_`sample' = 2*(1 - normal(t`var'_`sample'[1,1]))
		matrix  n`var'_`sample' = e(N)
		matrix           `var'tw_`sample' = [md`var'_`sample',p`var'_`sample']
		matrix  colnames `var'tw_`sample' =  mdtw`var'_`sample' ptw`var'_`sample'
	}
	matrix `var' = [.]
 	foreach sample of numlist 1(1)9 { 
		matrix `var' = [`var',`var'_`sample']
	}
	foreach sample in pool male female { 
		matrix `var' = [`var',`var'tw_`sample']
	}
	matrix `var' = `var'[1,2...]
}

// plot 
clear
foreach var in educ18 beh18 all18 educ18r beh18r all18r {
	svmat     `var', names(col)
	foreach sample of numlist 1(1)9 {
		gen t`var'_`sample' = md`var'_`sample' + c`var'_`sample'
	}
}

// labels
foreach var of varlist md* {
	tostring `var', gen(`var'_str) force format(%15.2fc)
	replace  `var'_str = "{&Delta} = " + `var'_str
}
foreach var of varlist n* {
	tostring `var', replace force format(%15.0fc)
	replace  `var' = "N = " + `var'
}
foreach var of varlist ptw* { 
	tostring `var', gen(`var'_str) force format(%15.2fc)
	replace  `var'_str = "{it:p} = " + `var'_str
}
gen h0ts   = "H{sub:0}: {&Delta}{sup:Singletons} = {&Delta}{sup:Twins}"
gen h0tsy  = .915
gen h0tspy = .895

gen h0tsyr  = .58
gen h0tspyr = .53

foreach num of numlist 1(1)85 {
	gen n_`num'  = `num'
	gen n_`num'p = n_`num' + .5
}


cd $output
// // by twinning
#delimit
twoway (bar     ceduc18_1  n_1, color(gs10) barwidth(.8)             )
	   (bar     teduc18_1  n_1, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter teduc18_1  n_1                                           , msymbol(none) mlabel(mdeduc18_1_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter teduc18_1  n_1 if   peduc18_1   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter teduc18_1  n_1 if   peduc18_1 >= 0.05 &  peduc18_1 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     cbeh18_1  n_5, color(gs10) barwidth(.8)             )
	   (bar     tbeh18_1  n_5, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tbeh18_1  n_5                                           , msymbol(none) mlabel(mdbeh18_1_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tbeh18_1  n_5 if   pbeh18_1   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tbeh18_1  n_5 if   pbeh18_1 >= 0.05  &  pbeh18_1 <= .10, msize(small)  mcolor(gs0) msymbol(square))

	  
	   (bar     call18_1  n_9, color(gs10) barwidth(.8)             )
	   (bar     tall18_1  n_9, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18_1  n_9                                           , msymbol(none) mlabel(mdall18_1_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tall18_1  n_9 if   pall18_1   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tall18_1  n_9 if   pall18_1 >= 0.05  &  pall18_1 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   

	   (bar     ceduc18_2  n_2, color(gs10) barwidth(.8)             yline(.55, lpattern(solid) lcolor(gs14)))
	   (bar     teduc18_2  n_2, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter teduc18_2  n_2                                           , msymbol(none) mlabel(mdeduc18_2_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter teduc18_2  n_2 if   peduc18_2   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter teduc18_2  n_2 if   peduc18_2 >= 0.05 &  peduc18_2 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     ceduc18_3  n_3, color(gs10) barwidth(.8)             )
	   (bar     teduc18_3  n_3, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter ceduc18_3  n_3                                           , msymbol(none) mlabel(mdeduc18_3_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ceduc18_3  n_3 if   peduc18_3   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter ceduc18_3  n_3 if   peduc18_3 >= 0.05 &  peduc18_3 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     cbeh18_2  n_6, color(gs10) barwidth(.8)             )
	   (bar     tbeh18_2  n_6, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tbeh18_2  n_6                                           , msymbol(none) mlabel(mdbeh18_2_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tbeh18_2  n_6 if   pbeh18_2   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tbeh18_2  n_6 if   pbeh18_2 >= 0.05  &  pbeh18_2 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     cbeh18_3  n_7, color(gs10) barwidth(.8)             )
	   (bar     tbeh18_3  n_7, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter cbeh18_3  n_7                                           , msymbol(none) mlabel(mdbeh18_3_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter cbeh18_3  n_7 if   pbeh18_3   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter cbeh18_3  n_7 if   pbeh18_3 >= 0.05  &  pbeh18_3 <= .10, msize(small)  mcolor(gs0) msymbol(square))

	  
	   (bar     call18_2  n_10, color(gs10) barwidth(.8)             )
	   (bar     tall18_2  n_10, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18_2  n_10                                           , msymbol(none) mlabel(mdall18_2_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tall18_2  n_10 if   pall18_2   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tall18_2  n_10 if   pall18_2 >= 0.05  &  pall18_2 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   
	   (bar     call18_3  n_11, color(gs10) barwidth(.8)             )
	   (bar     tall18_3  n_11, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18_3  n_11                                           , msymbol(none) mlabel(mdall18_3_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter call18_3  n_11 if   pall18_3   < 0.05                    , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter call18_3  n_11 if   pall18_3 >= 0.05  &  pall18_3 <= .10 , msize(small)  mcolor(gs0) msymbol(square))
	   
	   (scatter h0tspy  n_2p                                           , msymbol(none) mlabel(ptweduc18_pool_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter h0tsy   n_2p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter h0tspy  n_6p                                           , msymbol(none) mlabel(ptwbeh18_pool_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter h0tsy   n_6p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter h0tspy  n_10p                                           , msymbol(none) mlabel(ptwall18_pool_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter h0tsy   n_10p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   ,
	   
	   
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(small) position(6) region(lcolor(gs0)))
				 
		  ylabel(.55[.1].95, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Pooled" " " " " "' 2 `" "Singletons" " " " " "' 3 `" "Twins" " " " " "'
			     5 `" "Pooled" " " " " "' 6 `" "Singletons" " " " " "' 7 `" "Twins" " " " " "'
				 9 `" "Pooled" " " " " "' 10 `" "Singletons" " " " " "' 11 `" "Twins" " " " " "'	  
				 2 `" " " "{bf:Educational Index}" "' 6  `" " " "{bf:Behavioral Index}" "' 10  `" " " "{bf:Average of Indices}" "',  labsize(small) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export ncskills_bytwinning.eps, replace


#delimit
twoway (bar     ceduc18r_1  n_1, color(gs10) barwidth(.8)             yline(0, lpattern(solid) lcolor(gs14)))
	   (bar     teduc18r_1  n_1, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter teduc18r_1  n_1                                           , msymbol(none) mlabel(mdeduc18r_1_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter teduc18r_1  n_1 if   peduc18r_1   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter teduc18r_1  n_1 if   peduc18r_1 >= 0.05 &  peduc18r_1 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     cbeh18r_1  n_5, color(gs10) barwidth(.8)             )
	   (bar     tbeh18r_1  n_5, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tbeh18r_1  n_5                                           , msymbol(none) mlabel(mdbeh18r_1_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tbeh18r_1  n_5 if   pbeh18r_1   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tbeh18r_1  n_5 if   pbeh18r_1 >= 0.05  &  pbeh18r_1 <= .10, msize(small)  mcolor(gs0) msymbol(square))

	  
	   (bar     call18r_1  n_9, color(gs10) barwidth(.8)             )
	   (bar     tall18r_1  n_9, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18r_1  n_9                                           , msymbol(none) mlabel(mdall18r_1_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tall18r_1  n_9 if   pall18r_1   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tall18r_1  n_9 if   pall18r_1 >= 0.05  &  pall18r_1 <= .10, msize(small)  mcolor(gs0) msymbol(square))



	   (bar     ceduc18r_2  n_2, color(gs10) barwidth(.8)             yline(.6, lpattern(solid) lcolor(gs14)))
	   (bar     teduc18r_2  n_2, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter teduc18r_2  n_2                                           , msymbol(none) mlabel(mdeduc18r_2_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter teduc18r_2  n_2 if   peduc18r_2   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter teduc18r_2  n_2 if   peduc18r_2 >= 0.05 &  peduc18r_2 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     ceduc18r_3  n_3, color(gs10) barwidth(.8)             )
	   (bar     teduc18r_3  n_3, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter ceduc18r_3  n_3                                           , msymbol(none) mlabel(mdeduc18r_3_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ceduc18r_3  n_3 if   peduc18r_3   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter ceduc18r_3  n_3 if   peduc18r_3 >= 0.05 &  peduc18r_3 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     cbeh18r_2  n_6, color(gs10) barwidth(.8)             )
	   (bar     tbeh18r_2  n_6, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tbeh18r_2  n_6                                           , msymbol(none) mlabel(mdbeh18r_2_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tbeh18r_2  n_6 if   pbeh18r_2   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tbeh18r_2  n_6 if   pbeh18r_2 >= 0.05  &  pbeh18r_2 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   (bar     cbeh18r_3  n_7, color(gs10) barwidth(.8)             )
	   (bar     tbeh18r_3  n_7, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter cbeh18r_3  n_7                                           , msymbol(none) mlabel(mdbeh18r_3_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter cbeh18r_3  n_7 if   pbeh18r_3   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter cbeh18r_3  n_7 if   pbeh18r_3 >= 0.05  &  pbeh18r_3 <= .10, msize(small)  mcolor(gs0) msymbol(square))

	  
	   (bar     call18r_2  n_10, color(gs10) barwidth(.8)             )
	   (bar     tall18r_2  n_10, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18r_2  n_10                                           , msymbol(none) mlabel(mdall18r_2_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter tall18r_2  n_10 if   pall18r_2   < 0.05                   , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter tall18r_2  n_10 if   pall18r_2 >= 0.05  &  pall18r_2 <= .10, msize(small)  mcolor(gs0) msymbol(square))
	   
	   
	   (bar     call18r_3  n_11, color(gs10) barwidth(.8)             )
	   (bar     tall18r_3  n_11, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18r_3  n_11                                           , msymbol(none) mlabel(mdall18r_3_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter call18r_3  n_11 if   pall18r_3   < 0.05                    , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter call18r_3  n_11 if   pall18r_3 >= 0.05  &  pall18r_3 <= .10 , msize(small)  mcolor(gs0) msymbol(square))
	   
	   (scatter h0tspyr  n_2p                                           , msymbol(none) mlabel(ptweduc18r_pool_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter h0tsyr   n_2p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter h0tspyr  n_6p                                           , msymbol(none) mlabel(ptwbeh18r_pool_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter h0tsyr   n_6p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter h0tspyr  n_10p                                           , msymbol(none) mlabel(ptwall18r_pool_str) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter h0tsyr   n_10p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   ,
	   
	   
	   
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(small) position(6) region(lcolor(gs0)))
				 
		  ylabel(-.2[.2].6, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Pooled" " " " " "' 2 `" "Singletons" " " " " "' 3 `" "Twins" " " " " "'
			     5 `" "Pooled" " " " " "' 6 `" "Singletons" " " " " "' 7 `" "Twins" " " " " "'
				 9 `" "Pooled" " " " " "' 10 `" "Singletons" " " " " "' 11 `" "Twins" " " " " "'	  
				 2 `" " " "{bf:Educational Index}" "' 6  `" " " "{bf:Behavioral Index}" "' 10  `" " " "{bf:Average of Indices}" "',  labsize(small) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export ncskillsr_bytwinning.eps, replace

// by sex and twinning
#delimit
twoway (bar     call18_5  n_1, color(gs10) barwidth(.8)             yline(.55, lpattern(solid) lcolor(gs14)))
	   (bar     tall18_5  n_1, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18_5  n_1                                           , msymbol(none) mlabel(mdall18_5_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter tall18_5  n_1 if   pall18_5   < 0.05                   , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter tall18_5  n_1 if   pall18_5 >= 0.05  &  pall18_5 <= .10, msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   (bar     call18_6  n_2, color(gs10) barwidth(.8)             )
	   (bar     tall18_6  n_2, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18_6  n_2                                           , msymbol(none) mlabel(mdall18_6_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter call18_6  n_2 if   pall18_6   < 0.05                    , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter call18_6  n_2 if   pall18_6 >= 0.05  &  pall18_6 <= .10 , msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   
	   (bar     call18_8  n_4, color(gs10) barwidth(.8)             )
	   (bar     tall18_8  n_4, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18_8  n_4                                           , msymbol(none) mlabel(mdall18_8_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter tall18_8  n_4 if   pall18_8   < 0.05                   , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter tall18_8  n_4 if   pall18_8 >= 0.05  &  pall18_8 <= .10, msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   (bar     call18_9  n_5, color(gs10) barwidth(.8)             )
	   (bar     tall18_9  n_5, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18_9  n_5                                           , msymbol(none) mlabel(mdall18_9_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter call18_9  n_5 if   pall18_9   < 0.05                    , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter call18_9  n_5 if   pall18_9 >= 0.05  &  pall18_9 <= .10 , msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   (scatter h0tspy  n_1p                                           , msymbol(none) mlabel(ptwall18_male_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter h0tsy   n_1p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   
	   (scatter h0tspy  n_4p                                           , msymbol(none) mlabel(ptwall18_female_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter h0tsy   n_4p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))

	   ,
	   
	   
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(.55[.1].95, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Singletons" " " " " "' 2 `" "Twins" " " " " "'
			     4 `" "Singletons" " " " " "' 5 `" "Twins" " " " " "'
				1.5 `" " " "{bf:Male}" "' 4.5  `" " " "{bf:Female}" "',  labsize(medsmall) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export ncskills_bytwinningandsex.eps, replace


// by sex and twinning
#delimit
twoway (bar     call18r_5  n_1, color(gs10) barwidth(.8)             yline(.6, lpattern(solid) lcolor(gs14)))
	   (bar     tall18r_5  n_1, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18r_5  n_1                                           , msymbol(none) mlabel(mdall18r_5_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter call18r_5  n_1 if   pall18_5   < 0.05                   , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter call18r_5  n_1 if   pall18_5 >= 0.05  &  pall18_5 <= .10, msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   (bar     call18r_6  n_2, color(gs10) barwidth(.8)             )
	   (bar     tall18r_6  n_2, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18r_6  n_2                                           , msymbol(none) mlabel(mdall18r_6_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter call18r_6  n_2 if   pall18_6   < 0.05                    , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter call18r_6  n_2 if   pall18_6 >= 0.05  &  pall18_6 <= .10 , msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   
	   (bar     call18r_8  n_4, color(gs10) barwidth(.8)             )
	   (bar     tall18r_8  n_4, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter tall18r_8  n_4                                           , msymbol(none) mlabel(mdall18r_8_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter tall18r_8  n_4 if   pall18_8   < 0.05                   , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter tall18r_8  n_4 if   pall18_8 >= 0.05  &  pall18_8 <= .10, msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   (bar     call18r_9  n_5, color(gs10) barwidth(.8)             )
	   (bar     tall18r_9  n_5, fcolor(none) lcolor(gs0) barwidth(.8))
	   
	   (scatter call18r_9  n_5                                           , msymbol(none) mlabel(mdall18r_9_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter call18r_9  n_5 if   pall18_9   < 0.05                    , msize(medsmall)  mcolor(gs0) msymbol(circle))
	   (scatter call18r_9  n_5 if   pall18_9 >= 0.05  &  pall18_9 <= .10 , msize(medsmall)  mcolor(gs0) msymbol(square))
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   
	   (scatter h0tspyr  n_1p                                           , msymbol(none) mlabel(ptwall18r_male_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter h0tsyr   n_1p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   
	   (scatter h0tspyr  n_4p                                           , msymbol(none) mlabel(ptwall18r_female_str) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))
	   (scatter h0tsyr   n_4p                                           , msymbol(none) mlabel(h0ts) mlabposition(12) mlabcolor(gs0) mlabsize(medsmall))

	   ,
	   
	   
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(-.20[.20].60, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Singletons" " " " " "' 2 `" "Twins" " " " " "'
			     4 `" "Singletons" " " " " "' 5 `" "Twins" " " " " "'
				1.5 `" " " "{bf:Male}" "' 4.5  `" " " "{bf:Female}" "',  labsize(medsmall) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export ncskillsr_bytwinningandsex.eps, replace
